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It was recently shown by Jimbo, Miwa and Smirnov that the correlation 
. functions of a generalized XXZ chain associated with an inhomogeneous 

six-vertex model with disorder parameter a and with arbitrary inhomo- 
J ■ geneities on the horizontal lines factorize and can all be expressed in terms 

of only two functions p and go. Here we approach the description of the 
same correlation functions and, in particular, of the function co from a dif- 
ferent direction. We start from a novel multiple integral representation for 
the density matrix of a finite chain segment of length m in the presence of 
a disorder field a. We explicitly factorize the integrals for m = 2. Based on 
' this we present an alternative description of the function CO in terms of the 

solutions of certain linear and nonlinear integral equations. We then prove 
directly that the two definitions of co describe the same function. The def- 
inition in the work of Jimbo, Miwa and Smirnov was crucial for the proof 
■ of the factorization. The definition given here together with the known 

O description of p in terms of the solutions of nonlinear integral equations is 

useful for performing e.g. the Trotter limit in the finite temperature case, 
or for obtaining numerical results for the correlation functions at short dis- 
^ . tances. We also address the issue of the construction of an exponential 

Ctf \ form of the density matrix for finite a. 
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1 Introduction 

In recent years significant progress has been achieved in the understanding of the math- 
ematical structure of the correlation functions of the XXZ model and related integrable 
models. First of all the ground state correlation functions were studied. They are 
completely defined through the quantum-mechanical density matrix. An explicit ex- 
pression for the density matrix of a finite subchain of the infinite XXZ chain in the 
massive regime was first obtained by Jimbo, Miki, Miwa and Nakayashiki [17]. They 
expressed the elements of the density matrix in terms of multiple integrals. Subse- 
quently, extensions of their formulae to the massless regime and to a non-vanishing 
longitudinal magnetic field were obtained in [18,22]. 

Then it was realized that the multiple integrals can be factorized [10] and that, 
utilizing the so-called reduced Knizhnik-Zamolodchikov equation, the factorized in- 
tegrals can be written in a compact exponential form [6, 7]. The latter allows one to 
distinguish between an algebraic part and a physical part. The physical part is defined 
by a small number of transcendental functions, fixed by the one-point-correlators and 
by the two-point neighbour correlators which depend on the physical parameters like 
anisotropy, temperature, length of the chain, magnetic field, boundary conditions etc. 
The algebraic part is related to the representation theory of the symmetry algebra be- 
hind the model, namely the quantum group t/^sb) in case of the XXZ chain. 

In [8] it was observed that the formula for the correlation functions looks nicer 
if the XXZ chain is regularized by introducing an additional parameter, the disorder 
field a. With this new parameter it was possible to express the density matrix in terms 
of special fermionic annihilation operators b and c acting not on states of the spin 
chain, but on the space of (quasi-) local operators on these states. The annihilation 
operators appeared to be responsible for the algebraic part. The physical part was rep- 
resented by a transcendental function go determined by a single integral. In [9] the 
dual fermionic creation operators b*, c* and a bosonic creation operator t* were con- 
structed. These operators together generate a special basis of the space of quasi-local 
operators. Since b*, c* are Fermi operators, Wick's theorem applies and expectation 
values of products of b*, c* and t* in an appropriately defined vacuum state can be 
written as determinants, very much as in the case of free fermions. 

Thermodynamic properties of integrable lattice models can be studied within the 
Suzuki-Trotter formalism by considering an auxiliary lattice with staggering in the so- 
called Trotter direction [24]. The temperature appears as a result of a special limit 
when the extension of the lattice in Trotter direction becomes infinite. Physical quanti- 
ties are expressed in an efficient way through the solution to certain non-linear integral 
equations [23]. A detailed discussion of this issue and further references can, for in- 
stance, be found in the book [12]. 

In papers [14, 15] the Suzuki-Trotter formalism was used in order to generalize 
the multiple integrals to finite temperature. Then their factorization was probed for 
several examples of correlation functions, first for the XXX chain [3] and later for the 
XXZ chain [2,4]. Also a conjecture was formulated stating that the above mentioned 
exponential form is valid with the same fermionic operators (at least as long as they 
act on spin reversal invariant products of local operators) as for the ground state and 
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two functions go, go' obtained from an oc-dependent function in the limit a — > 0. 

Unfortunately, the formulae of [2, 4] worked only in this limit. The generalization 
to generic a stayed obscure. One of the purposes of the present work is to add to 
the clarification of this point, starting from a proper multiple integral representation 
with disorder parameter a. Here, as we had to learn [20], the crucial point is that the 
'Cauchy extraction trick', invented in [16] and described in more detail in [21], can be 
applied in the finite temperature case and also in the more general situation of a finite 
lattice with inhomogeneities in Trotter direction. 

Important new insight came from a recent paper [19] by Jimbo, Miwa and Smirnov, 
where they suggested a purely algebraic approach to the problem of calculating the 
static correlation functions of the XXZ model. The key idea of [19] is to evaluate a 
linear functional related to the partition function within the fermionic basis constructed 
in [9]. The authors of [19] work with a finite lattice, inhomogeneous in Trotter direc- 
tion. In this situation they suggest a new and surprising construction of the function go 
depending on a magnetic field and on the disorder parameter a. 

In the present paper we discuss the relation of the work by Jimbo, Miwa and 
Smirnov to the approach using non-linear integral equations which at the moment 
seems more appropriate e.g. for taking the Trotter limit (which was omitted in [19]). 
In particular, we present an alternative description of the function go starting from 
the multiple integral and using the explicit factorization of the density matrix for two 
neighbouring lattice sites. We then give a direct proof that our expression, though 
looking rather different than that in [19], in fact describes the same function. 

An inhomogeneous lattice in Trotter direction is very general and leaves many 
different options for the realization of physical correlation functions. Here we shall 
concentrate on two of them, the correlation functions of the infinite XXZ chain at 
finite temperature and magnetic field (temperature case) and the ground state correla- 
tion functions of a finite chain with twisted periodic boundary conditions (finite length 
case). Both cases can be treated to a very large extend simultaneously. They are only 
distinct in that a different distribution of inhomogeneity parameters is required and in 
that for the finite temperature case the Trotter limit has to be performed. Note that 
instead of the XXZ Hamiltonian we could consider combinations of conserved quan- 
tities obtained from the transfer matrix of the six- vertex model within the formalism 
of non-linear integral equations. For the bulk thermodynamic properties this issue was 
recently studied in [25]. 

The paper is organized as follows. In the next section we define our basic objects 
and recall some of their properties. In the third section we show the multiple integral 
formula for the elements of the (a-twisted) density matrix for a sub-chain of length m. 
In section four we consider the simplest case m = 1. The fifth section is devoted to 
applying the factorization technique to the double integrals for m = 2. In section 6 we 
introduce the function go. We discuss its properties and the relation to its realization 
by Jimbo, Miwa and Smirnov. The content of section 7 is some preliminary work on 
the construction of an operator t, dual to the creation operator t*, which should appear 
in the construction of an exponential form for finite temperature and finite disorder 
parameter. In the appendices we provide a derivation of the multiple integral formulae, 
we discuss the limit a — > 0, and we compare with the results of the papers [2,4]. 
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2 Density matrix and correlation functions 



The XXZ quantum spin chain is defined by the Hamiltonian 

H N (K) = J £ (o^O* + 0^ ^ • A!0 } ,0) - 1)) , 



(i) 



: el a 3 " = i(el — e + ), o z = e 



+ 



written here in terms of the Pauli matrices o x ■ 
(where the are the elements of the gl(2) standard basis). The two real parameters J 
and A control the ground state phase diagram of the model. For simplicity of notation 
we shall restrict ourselves in the following to the critical phase / > 0, |A| < 1. Note, 
however, that the results of this work can be easily extended to the off-critical anti- 
ferromagnetic phase A > 1 . We shall also assume without further mentioning that the 
number of lattice sites N is even. 

To fully specify H^(k) we have to define the boundary conditions. We shall con- 
sider twisted periodic boundary conditions, when we are dealing with the ground state 
of the finite chain. Then Hn(k) depends on an additional parameter K through 




q -K<? I e N% e N- \ q K a z 
\ e N+ e N ZJ 



(2) 



Here q is related to A as A = (q + q~ l )/2. For the finite temperature case we shall 
assume periodic boundary conditions for the Hamiltonian. Nevertheless the same pa- 
rameter K will appear in that case as a twist parameter of the quantum transfer matrix, 
having then a rather different physical meaning as an external magnetic field coupling 
to the spins by a Zeeman term. We shall elaborate on this below. 

The integrable structure behind the Hamiltonian (1) is generated by the trigono- 
metric /^-matrix of the six-vertex model [1], 

/l 0\ 

b{X) c{\) 

c(k) b{X) 

\0 1/ 



R(X) = 



(3) 



b(X) 



sh(X) 



c{X) 



sh(r|) 



(4) 



sh(X + ri)' v ; sh(X + r|)' 

acting on C 2 (g> C 2 . As presented here it satisfies the Yang-Baxter equation in addi- 
tive form. To facilitate the comparison with [9, 19], where the multiplicative form 
was preferred, we set q = e 11 and £ = e\ Then for arbitrary complex inhomogeneity 
parameters Py, j = 1, . . . ,N, the definition 

T a (Q = RajfiX ~ Pat) . ../?«,! (I - Pi ) (5) 

of the monodromy matrix makes sense, where, as usual, the indices l,...,N refer to 
the spin chain while a refers to an additional site defining the so-called auxiliary space. 
We also set T a (C,,K) = T a (Qq KG " a and introduce the twisted transfer matrix 



f(C,K)=tr B (r a (S,K)) 



(6) 
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In [19] a six vertex-model with N horizontal rows and an arbitrary distribution of 
the inhomogeneities ij = e^ on these rows was considered. Here we would like to 
point out that two specific distributions are of particular interest in physical applica- 
tions. Moreover, in both cases the special functions that enter the representations of the 
transfer matrix eigenvalues and correlation functions have nice descriptions in terms 
of solutions of linear and non-linear integral equations. 

The first case relates to the ground state of the Hamiltonian (1). We call it the finite 
length case. In this case we choose 

Py = Tl/2, j=l,...,N. (7) 

Then 

H N (K) =2/sh(ri)a x ln(r 1 (^i,K)?(C,K))| x=ii/2 , (8) 

with twisted boundary conditions (2) if we identify A = ch(r|). The critical regime 
|A| < 1 corresponds to purely imaginary r\ = ij, j <G [0,7c). In this case the physical 
twist angle or flux <1> G [0,2tc) is <I> = — icy, whence K should be real. If we stick to the 
vertex model picture of [19], then ?(£,k) is the vertical or column-to-column transfer 
matrix in this case. 

The second case is determined by an alternating choice 




P; = So S > J = h-,N/2, (9) 



of inhomogeneity parameters. This case will be called the finite temperature case as it 
relates to the quantum transfer matrix, whose monodromy matrix is 

T QTM {Q = 

R a , N & ~ P/A0^-i,«(-P/Ar - X) . . .R a , 2 (l - P/A^X(-P/iV - 1) . (10) 

Here the superscript 'ti ' indicates transposition with respect to the first space. In fact, 
setting Y = J7 l^y-i an d using the crossing symmetry 

a y j R aJ (X-n)a y j = b(X-n)R t l a (-X) (11) 
of the /^-matrix we find that 

N/2 1 

The quantum transfer matrix is by definition 

t^ M (M=tr a (l^ M {M), (13) 

where T?™ = T?™ (Q q ^« . 

Again, within the vertex model picture, t@™(C, ,k), or t(C,,K) with the choice (9) 
of inhomogeneity parameter, corresponds to the vertical transfer matrix. There is an 
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important difference, though, that has been explained at several occasions [14,23]. In 
the finite length case the Hamiltonian can be derived from the vertical transfer ma- 
trix. In particular, the vertical transfer matrix and the Hamiltonian (1) have the same 
eigenstates. In the finite temperature case, on the other hand, with a lattice which is 
homogeneous in horizontal direction, say, the Hamiltonian is related to the horizontal 
transfer matrix with purely periodic boundary conditions. It is then also periodic and 
will be denoted Hl(0), where L is the horizontal extension of the lattice. In this case 
the vertical transfer matrix eigenstates are different from those of the Hamiltonian. In 
particular, the eigenstate with the largest modulus determines the state of thermody- 
namic equilibrium in the thermodynamic limit, i.e. the free energy of the XXZ chain 
and all its static correlation functions [14]. Also the physical interpretation of the pa- 
rameter K is rather different in this case. It corresponds to a magnetic field coupling to 
the spin chain through a Zeeman term (see e.g. [14]). 

Using a lattice of finite extension L in horizontal direction we can express the 
partition function of the homogeneous XXZ chain of length L as 

Z L = tn LfS -P*dO)+hS M /T = i imt r li ... jJV ( f e™(l,V(2T 1 r))) L . (14) 

Here T is the temperature and h is a longitudinal magnetic field. P must be chosen as 
P = 2/sh(r|)/r. Furthermore 

^ = jt°5 < 15 > 

is the conserved z-component of the total spin. Equation (14) becomes efficient in the 
thermodynamic limit L — > °°, since then a single eigenvalue A^ rM (l,K) of t^™{\,K) 
of largest modulus dominates the large-L asymptotics of Zi in the Trotter limit N — > °°. 
We shall refer to this eigenvalue as the dominant one. 

We would like to remark that in our understanding the quantum transfer matrix 
is, in general, more appropriate for studying integrable spin models on the infinite 
lattice than the usual transfer matrix. In general there is no crossing symmetry, and the 
quantum transfer matrix and the usual transfer matrix are not related by a similarity 
transformation like in (12). Also within the quantum transfer matrix formulation the 
density matrix directly takes its 'natural form' in terms of monodromy matrix elements 
(see below). No solution of a quantum inverse problem as in [22] is required. In our 
particular case we do have the crossing symmetry, and the quantum transfer matrix 
and the usual transfer matrix with staggered inhomogeneities (9) give an equivalent 
description of the density matrix (see below). Still, the largest eigenvalue of t(C,,K) 
with the distribution (9) of inhomogeneities diverges in the Trotter limit as can be seen 
from (12). 

Let us come back to the situation of arbitrarily distributed inhomogeneity param- 
eters P ; . Following [19] we shall assume that for a certain spectral parameter i^o and 
any K € C the transfer matrix f(^o> K ) has a unique eigenvector |k) with eigenvalue 
A(^o,k) of largest modulus. This is certainly true for the two special cases considered 
above. In the finite length case i^o = l 1 ^ 2 , while £o = 1 in the finite temperature case. 
We fix a set of 'vertical inhomogeneity parameters' Vi, . . . ,V m and set \j = eV Then 
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we can define the object of our main interest, the density matrix with matrix elements 



D We ,.. em (^,...,^|K,a)- (K + a|n7 =1 r(^,K)|K) ' (16) 

which is in fact an inhomogeneous and 'oc-twisted' version of the usual density matrix. 

In the finite length case (7) with twist angle <I> the expectation value in the ground 
state |<J>) of any operator acting non-trivially only on the first m lattice sites 

is [11] ' ' 

^ii^ = So v l^^-ADN&U • • • ~ 3>/Y,0C)*[l,m] } ■ (17) 

In the finite temperature case (9) we use that the right hand side of (16) stays form in- 
variant under the transformation (12) which replaces all objects relating to the ordinary 
transfer matrix with the corresponding objects relating to the quantum transfer matrix. 
Hence, from [15], 

,„ . tr 1 ,.,. x {e-^(o)^ |1 . 1 /^ [1 , m] } 
{X[i,m])T,h = hm — 

= lim lim limtn m {D N (^, . . . ,^ m \h/(2^),a)X [Lm] } . (18) 

The density matrix (16) allows for reduction from the left and from the right ex- 
pressed by 

tri {Dn^ , . . . , £ m |K, a)q a(f i } = p(£i )D N (^ 2 , ■ ■ ■ , ^»|k, a) , (19a) 
tr m {D JV (^i,...,^ m |K,a)} =£) JV (^i,...,^ m _i|K,a), (19b) 

where 

The function p plays an important role in [19]. As we shall see below it is also im- 
portant for the formulation of a multiple integral formula for the density matrix and 
is the only non-trivial one-point function for finite a. In the temperature case with 
K = h/(2r\T) we have 

p(l) = l+m(7»2r|a + 0(a 2 ), (21) 

where m(T,h) is the magnetization. 

In the temperature case as well as in the finite length case and in certain inhomo- 
geneous generalizations of both cases the function p can be expressed in terms of an 
integral over certain auxiliary functions (see e.g. [11, 14]), 



1 + a( i u,K + oc) 



(22) 



'e 27ti [ l + a^K) 

Here e(X) is the 'bare energy' 

e{X) = cth(X) - cth(A, + r|) (23) 
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dm 



Re 



Figure 1 : The canonical contour S surrounds the real axis in counterclockwise manner 
inside the strip — \ < ImX < \. 



and o(X, k) is the solution of a non-linear integral equation with integration kernel 

AT(A,) = cth(A,-Ti)-cth(A, + 'n). (24) 
In the finite length case this equation reads 

ln(c(A,,K)) = 



N 



(N - 2k)t\ + £ In 



r sh(x-p 7 ) 




sh(X-p y - + ri)_ 


-L 



d^ 



'e 27ii 



7 ff(A,-/i)ln(l + o(/i,K)). (25) 



Equations (22) and (25) are still valid if the Py are not precisely those of equation (7), 
but are close to r|/2 with Im (3 ; - = y/2. The contour of integration to be used in (22) and 
(25) is shown in figure 1. In the temperature case the non-linear integral equation has 
a similar structure, but the driving term is different. Suppose that for j = l,...,N/2 
the P27-1 are close to r\, whereas the $2j are close to 0. Then 



ln(a(^,K)) = -2kt| 

^ 2 f sh(A,-p 2 j)sh(A,-p 2 ;- 1 +2Ti) 

^ L 11 sh(X-p 27 -+Tl)sh(X-p 2y _ 1 +Tl)J 



| e ^(X-/i)ln(l + aCu,K)). (26) 



We presented both equations (25) and (26) in inhomogeneous form, since we shall 
need this later, when comparing with [19]. Note, however, that the homogeneous limit 
is trivial in both cases and that, moreover, the Trotter limit can be performed in (26). 
Then 



ln(a(X,,K)) = (N-2K)r|+Mn 
in the finite length case and 
ln(o(A,,K)) = -2kt| - 



sh(X,-ri/2) 



^|jU(A,-/i)Ml + aOi,K)). (27) 



sh(X + r|/2)_ 
27sh(ri)e(X) f d^i 



f Q ^K{X-n)\n{\ + a{n,K)) (28) 
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in the temperature case and in the Trotter limit. Equations (27) and (28) are what 
we call the o-form of the non-linear integral equation. There is another so-called bb- 
form [11, 23] which is more convenient for an accurate calculation of the numerical 
values of the functions. 



3 The multiple integral representation of the density matrix 

In appendix A we derive the following multiple integral representation for the elements 
of the density matrix. 

' ' r p r i r m r 

D^\:l:& u ...£ m \K,a) = n / dm(Xj)F+(Xj) n / dm (h)F f -(Xj 

del M : m \ G(Xj,V k )] 



rii<y<*<msh(^ - X k - n) sh(v fc - v 7 ) 

where we have used the notation 



(29) 



dm ^ = 9,infnnLi dm(X) = a(X,K)dm(X), (30) 
27lip(y(l + a(A,K)) 



tj-1 m 

FHX)=Y\M^-Vk) J! sh(A.-v*=FTi), £j 

1 k=\ k=£ J + \ { £ m-j+l J = P+l,---,m 



4 j=l,...,p 



with 8+ the jth plus in the sequence (£j)J =i , £j the jth minus sign in the sequence 
( £ 'j)J=i and P tne number of plus signs in (£j)J =l . The function G is new here. It is 
defined as the solution of the linear integral equation 

G(X,v) =^ a cth(X-v-r|)-p(^)cth(X-v)+ f dm(v)K a (X-v)G(v,v), (31) 
where ^ = e v , and the kernel 

K a {X)=q- a cth(X-\\)-q a cth{X + r [ ) (32) 

is a deformed version of (24). 

Equation (29) is a generalization to finite a of the multiple integral formulae first 
derived in [11, 15]. To simplify the notation we shall sometimes suppress the depen- 
dence of the density matrix elements on K and a. 



4 The case m = 1 

For m = 1 there are only two non-vanishing density matrix elements. They are related 
to the function p by the reduction relations (19) which imply that 




1 -q 



-a 



q a _ q -a I q a 




(33) 
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When we insert equation (29) for m = 1 here, we do not obtain an independent equa- 
tion, but rather an interesting identity for p (recall that p appears in the measure), 

p© = - (q a - q~ a ) f dm(^G^v) . (34) 

Je 

It allows us to calculate the asymptotic behaviour of the function G, 

lim G(X,v)=0. (35) 

5 Factorization of the density matrix for m = 2 

The factorization of the multiple integrals for the ground state density matrix was dis- 
covered in [10]. In that case the integrand consists of explicit functions whose analytic 
properties were used in the calculation. In the finite temperature case a different fac- 
torization technique had to be invented. As was demonstrated in [3] the linear integral 
equation for the function G, appropriately used under the multiple integral, can be 
viewed as the source of the factorization, at least for the special case of the isotropic 
chain at a = 0. For the XXZ chain outside the isotropic point and without the disorder 
parameter a, however, that trick does not work anymore. Here we shall see that a finite 
a allows us to perform the factorization of the density matrix in much the same way 
as in [3]. 

Let us consider m = 2 in (29). There are six non-vanishing matrix elements in this 
case, one for p = 0, four for p = 1 and one for p = 2. We shall concentrate on the case 
p = 1, since the matrix elements for p = or 2 can be obtained from those for p = 1 
by means of the reduction relation (19). After substituting wj = e 2 ^ and t,j = e v J , 
j = 1,2, the corresponding integrals are all of the form 

J= g2 1 g2 f dm(Mi) f dm(fi2)det[G(fij,\ k )]r(wuW2), (36) 
q 2 -Si^e Je 

where 

p(wi , W2) 

r(w\,w 2 ) = s — , p(wi,w 2 ) =c w\W2 + c\W\ +C2W2 + C3. (37) 

w\ — q L W2 

The coefficients cj are different for the four different matrix elements. They are listed 
in table 1 . 
Inserting 

ds( "» = ^)- dmM (38) 

into (36) and taking into account that p(e^) is analytic and non-zero inside C we obtain 

5 (5i-s)=-/d mW de,( G ^r 2 \ ) 

- J dm (jj 1 ) dm (fi 2 ) det [G (jij , V* ) ] r ( w 1 , w 2 ) , (39) 
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e, e 2 

El &2 




Cl 


C 2 


C3 


+ - 

+ - 

- + 

- + 

+ - 

- + 

- + 
+ - 


1 

q 2 
qki/^2 


-qt>iZ>2 

-q-%^2 


~q 2 % 

-qt>iZ>2 

-q%Z,2 


q 2 m 
mi 

q^2 
q^l 



Table 1: The coefficients of the polynomial p. 



where w = e 2 ^. Here the first integral is already factorized. Under the second integral 
the integration measures now appear symmetrically. Hence, we may replace r(w\,w 2 ) 
by (r(wi,w 2 ) -r(w 2 ,wi))/2. 

Following [3] we want to use the integral equation (31) under the second integral 
in (39). This is possible if rational functions F(wi,w 2 ) and g(w) exist, such that 

r(wi,w 2 )-r(w 2 ,wi) = F(wi,w 2 ) + g(wi)K a (jj 1 -/j 2 ) -g(w 2 )K a (jj 2 -jui) , (40) 

and the antisymmetric function F(w\,w 2 ) is a sum of factorized functions in w\ and 
w 2 . Then F considered as a function of w\ cannot have poles whose position depends 
on w 2 . In particular, the residue at w\ = q 2 w 2 must vanish. Using this in (40) with the 
explicit forms of r and K a inserted we obtain a difference equation for g, 

(2 \ — i / \ p(q 2 w,w) 

g{qw)y -g(w)y= ^ . (41) 

Here y = q a . Clearly this equation has a solution of the form 

g_ 

g(w) =g+w + g -\ . (42) 

w 

The coefficients are easily obtained by substituting the latter expression into (41), 

coy c 3 y _ (ci+^ 2 c 2 )y 

g+ 2{q 2 -yiy 8 - 2(1-^2)' 80 2(l-y2) " ^ 

Substituting g back into (40) we obtain F(w\,w 2 ) = f{w\) —f(w 2 ), where 

f(w) = {y-y- l )(g + w- g —). (44) 

Consequentially 

r(wi,w 2 ) = f{w\) + g{w\)K a (fi\ —fi 2 ) + symmetric function. (45) 



With this we can factorize the second integral in (39) by means of the integral 
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equation (31), 

j^m{n\ ) f dmfa) det [G(jij,V k )] r(w\ , w 2 ) 

= , _ -l ) det ( v - S-9- (vi ) £+<P+ (v 2 ) - g-9- 

\ 9o(vi) 9o(v 2 ) / 

+ fdm( M )6t t ( , G( ^ Vl) 1A , G( ^ V2) 1 V (46) 
Je W [gWHfavw- 1 ) g(w)H(^v 2 ;y- 1 )) 

where 

cp,-(v) = / dmMwJGfav) , j = +,0,- , (47a) 
Je 

H( M ,v;y- l ) = p(Z > )cth( M -v)-y- 1 C th( M -v-r\). (47b) 

Finally we substitute (46) into (39) and further simplify the resulting expression using 
the identities 

gWHfaviy- 1 ) = gfflHfaviy) - cth(ju - v - t,) 

- (po(v) (f(w) - fg 2 )) + {m 2 ) - f{q 2 Z, 2 )) , (48a) 



r ^- PJ W lct ^- v - n) - p± W 1 - (48b) 



Then 



(y- 1 -p(^i))(y-pfe))/(^)-(y- 1 -pfe))(y-pfti))/(^) ™ 
(g-g)Cy-^ 1 ) 2 ' 



(50) 



where 

^1,^2) = | e dmCu)G( / i,V 2 )(^cthOu-v 1 -Tl)-p(^ 1 )cthCu-V 1 )). 

Equation (49) determines the four density matrix elements for p = 1 in factorized 
form. Note that the matrix elements depend on only two transcendental functions p 
and m. The remaining two non-zero density matrix elements for m = 2 follow from 
(49) by means of the reduction relations (19), 

D++(^ 2 ) = P ^ l) _~_\ l -Dlz&u^), (51a) 

DZZ(^ 2 ) = y -0^--D-_t(^2). (51b) 

We shall give a fully explicit matrix representation of the factorized density matrix for 
m = 2 below, after we have introduced the function oo. 
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6 The function 0) 

In the recent work [19] is was shown that the correlation functions defined by the 
inhomogeneous and oc-twisted density matrix (16) factorize and can all be expressed 
in terms of only two transcendental functions, the function p entering the reduction 
relations (19) and another function go which in [19] was defined as the expectation 
value of a product of two creation operators and was represented by a determinant 
formula. The approach of [19] is slightly different from ours here in that the lattice 
used in [19] is homogeneous in 'horizontal direction' (all the ^s in (16) are taken to 
be 1 from the outset). For the ground state both cases lead to the same function go 
(see section 5.3 and 5.4 of [9]). In particular, in the inhomogeneous case, following 
sections 5.1 and 5.3 of [9], we have* 

co(^^ 2 ) = -<cf 1)2] fe,a)bf 1)2] (^,a-l)(l)>. (52) 

Replacing the vacuum expectation value by the expectation value calculated with 
the density matrix (16) we take (52) as our definition of the function ox In our case go 
depends on two twist parameters K and a. We indicate this by writing go(^i,^2|k,oc). 
The construction of the operators b*j 2 ] and c*j 2 j is explained in [9]. For the product 
needed in (52) we find the explicit expression 

^ a c^ 2] fe,a)b[ lj2] (^,a-l)(l) = 

^x-l^-l q 1 -^- 1 ^q a -q- a 



where £ = ^1/^2- Inserting this into (52) and calculating the average with the factor- 
ized two-site density matrix of the previous section we obtain 

C0(^ 1 ^ 2 |K,a)=2^ 1 ^ 2 )-A ¥ © + 2(p(^ 1 )-pfe)) ¥ ©. (54) 
Here we adopted the notation from [9], 

and A is the difference operator whose action on a function / is defined by A/(^) = 

/(<£)- /(rt). 



^More precisely this function was denoted ((Oq — /£,2, °0 in [9]. 
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The remaining part of this section is devoted to the exploration of the properties 
of CO. First of all we substitute go back into the equation for the two-site density matrix, 
which can then be expressed entirely in terms of go and a function 

ch(om) - p(£) 

*«'* a >- sMon) (56) 

which is sometimes more convenient than the function p itself. We obtain 

Div(^i,^2|K,a) = -I 2 ®h 



4(^«-i-^-«)V %-%~ x 



|^-/2®O Z -^y-O Z 0(f+r 1 C + ®(J-+^®O + J 
1 ^- a - 1 C0 12 - ^ +1 0>2i + <?i<?2(q a -q- a ) 



(^- q ^ h®G z - <1+ 2 g z ®g z + £a + ®a~ + £~ 1 a~(g)a + ^ 

^ (cpi a z ® 7 2 + cp 2 7 2 (8) g z ) - 4 ^ == | rT j ((pi - 92) (a + 8g'-o~®c + ), (57) 



where we introduced the abbreviations GOyyt = Go(^y,^|K,a) and q> ; - = q>(^|K,oc). 

For the limit a — > the properties of the functions (p and go with respect to negating 
K and a are important. They follow from the fact that the /^-matrix is invariant under 
spin reversal, 

R(l) = (<?®<?)R(X){& c ®<f). (58) 
Introducing the spin reversal operator J = a\ . . . we conclude with (58) that 

T a ^,-K)=a x a JT a ^,K)J(f a . (59) 

It follows that t(C„ -k) = Jt(^,K)J. Hence, 

/|k) = |-k), (60a) 

A(C,K)=A(C,-K). (60b) 

The latter two equations used in the definition (16) of the a-twisted density matrix 
imply that 

D N g u ...,U-^-Oi) = (o x r m D N ^ u ...,U^)(o x f m . (61) 
From (20), (60b) we obtain the relation 

q>(q-K,-a) = -<p(C|K,a). (62) 
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Equation (61) together with (52)-(54) and the expressions for the density matrix ele- 
ments of the previous section implies that 

(o(£i,§2|K,a) = a>($2,5i|-K,-a). (63) 

Our next step is to verify that the function go given by the formula (54) satisfies 
a property called the 'normalization condition' by the authors of [19] (see equation 
(6.10) there). So we come back to the case of finite Trotter number N with arbitrary 
inhomogeneity parameters P 7 , j = 1, . . . ,N as it is written in (5). We shall also use 
multiplicative parameters Xj = e$' . 

We consider the normalization condition in the following form 

(co(C,^|K,a)+D^-V(CA))| c=T .+ 

+ p(x;)(©(C^|K,a)+D c D^V(CA))| c=9 -i t . =0, (64) 

j = 1, . ..N, which can be obtained from the integral in (6.10) of [19] by taking the 
residues and using the TQ-relation (4.2) of that paper. Also let us recall the definition 

D^g(Q=g(qQ+g(q- 1 Q-2p(Qg(Q. (65) 

Actually, (6.10) of [19] comprises one more equation related to the residue at ^ 2 = 0. 
This case needs separate treatment and will be discussed below. 

First we use the following difference equation for the function *P defined by (50), 

+ pfe)cth(vi -v 2 )-^ a cth(vi -v 2 -r|) 

-q- a {p(^)p(q-%)-l) f dmCu)G(/i,V 2 )cth( / i-V 1 +Tl), (66) 

J c 

where a = l/o by definition. This equation is the result of an analytical continuation 
defined for l*^ -1 ^,^) through an appropriate deformation of the integration contour 
in (50). Some simplifications occur in the limit Vi — > [3/ or equivalently ^ — > X/, 
namely, since a((3 ; -,K) = a($j — T|,k) = or 6(13/, k) = a(f3 ; - -T|,k) = oo, the first two 
terms in the right hand side of (66) do not contribute. Then we have 

n(r w -u \ = 6~(g~ 1 Tj;K + a)j2 + (T y -;K) 6-(T y -;K + q)j2 + (g~ 1 T J -;K) = 
PUiPW j) Q -( T . ;K+a ) e +^-i T . ;K )- e -( ? -i T . ;K + a)e+ ( T . ;K ) 

with the Q-f unctions Q ± defined in [19]. This means that also the last term in the right 
hand side of (66) does not contribute. Hence, we obtain 

^(x/,^ 2 ) + p(x^- a ^- 1 x i ^ 2 )=pfe)cth(P y -v 2 )-^ a cth(P 7 -v 2 -Ti). (67) 

Note that the right hand side is, up to the sign, equal to the driving term in the integral 
equation (31) for G. 
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If we take the formula (54) and use (67) then, after some algebra, we obtain 

co(t;, §2|k, a) + p (t;)w(?~S-> ^|k, a) = 

- (A^(^ 2 )) |^ - p(x 7 ) (^(C/fe)) |^ 

+ 2(p(x y ) + pfe)) ¥ (x i -A 2 )-2(l + p(x y )pfe)) ¥ (^ 1 x i -A 2 ). (68) 

Now we need to check that this equation is equivalent to (64). To this end we should 
verify the following equality 

(D^V(C/5))|^ + P(Ty)(D c ^V«/«)|^ = 

(a ;¥ (CA 2 ))|^ + p(x,)(a c¥ (^ 2 ))| c=? _ 1Xj . 

-2(p(x i ) + pfe)) ¥ (x i -A 2 )+2(l + p(x i -)pfe)) ¥ (,?- 1 x y A 2 ). (69) 

Using the definition (65) we come after a little algebra to the following expression for 
an arbitrary function g(Q 

d& = +4(1 - P (Q) (i - P ft)) g (w 

-2(p(Q + p©) (g(^A) + g(?- 1 CA)-2g(CA)). (70) 



Now take 



(D^g(^))\^+p(Tj)(D^ g (^)) 



(A\ g(C/£)) 1^ + (A* |^ - 2(p(xy) + pft)) (Ac *(C/«) 1^ 

+ 2(l+p(x y )p©)(A^(CA))|^_ lT .. (71) 

If we substitute = A^ 1 \|/(^/^) and take £, = ^ 2 , then we immediately arrive at 

the equality (69). 

As was mentioned above, there is one more case to be considered, corresponding 
to the contour To, i.e. to the residue at ^ 2 = in equation (6.10) of [19] which has to 
vanish. Its vanishing follows from 



2q- K 



q K + q K 



Pfe)- q- a + (q a -q- a ) [ dm(p)G(^V 2 ) 

Je 



0. (72) 



Here we have used (50), (54), (55), (70) as well as the fact that limv-»-<»p(£) = 
(q a+K + q~ a ~ K ) I (q K + g~ K ) in the first equation and the identity (34) in the second 
equation. 

The normalization condition just shown to be satisfied by our function go defined in 
(52) is the main ingredient in our proof that go is in fact the same function as introduced 
in equation (7.2) of [19]. Let us consider go as a function of As was shown in [19] 
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Figure 2: Four cases to be considered for the analytic continuation of ¥(£ 1,^2) with 
respect to Vi. Here C is the canonical contour of figure 1. 



the function p(^i) depends only on £j. The same is then true for ^(^1,^2) from (50). 
Using (54) we conclude that ^~ a Go(^i,^ 2 |K,a) is a function of £j. From its definition 
(52) and from (16), (53) we see that GO is rational in ^ of the form P(£j) /Q(£i)> where 
P and 2 are polynomials. Clearly both of them are at most of degree TV + 2. The zeros 
of Q are the N zeros of the transfer matrix eigenvalue A(^i ,k) plus two zeros at q ±2 i^ 
stemming from the two simple poles of ^ _a cfj 2 , (^Ot)!}^ 2 ,(^i,a — 1)(1). Comparing 
now with the definition (7.2) of [19] we see that the functions there has precisely the 
same structure. It is rational of the form P(£j)/<2(£j) with two polynomials P, Q at 
most of degree N + 2. Q and Q have the same zeros. We may therefore assume that 
they are identical. In order to show that P and P also agree we have to provide N + 3 
relations. N + 1 of them are given by the normalization condition above. Another two 
come from the residues at the two trivial poles. 

Since they are outside the canonical contour, we have to consider again the analytic 
continuation of the integral (50) defining *P with respect to £1. There are four regions 
depending on the location of Vi relative to the contour (see figure 2). Using (50) we 
obtain 

V$i,Z t2 ) = I dm( / i)GCu,v 2 )(^cth( / i-v 1 -Ti)-p^ 1 )cthCu-v 1 )) 

r g(vi,v 2 ) _ 

case (I) 

1 + o(Vi,k) 

case (II) 

G(v„v 2 ) ," C (v, +1) v 2 ) ^ 



l + o(v b K) (1 + o(vi+ti,k))p(^i) 

G(Vi,V 2 ) 

k 1 + o(vi,k) 



case (IV) . 



18 



Then, e.g. be means of the integral equation (31) 



^f-*a(% 2 q)d(% 2 ) 



(l + a(v 2 +ri,K))(l + o(v 2 ,K)) A(^,K)Afe,K) ' 



( 1 + a(v 2 , k) ) ( 1 + a(v 2 - ri , k) ) Afe , k) A{% 2 q- 1 , k) 



(74a) 



(74b) 



where a and d are the vacuum expectation values of the diagonal elements of T(Q. 
Since the ratios on the right hand side are invariant under changing the normalization 
of the /^-matrix, we can directly compare the residues obtained from (54), (74) with 
those obtained from equation (7.2) of [19]. We find agreement, which completes the 
proof. 

What if we consider (54) as the definition of CO? Then, in addition, we have to 
show that there are no poles other than the two trivial ones and those at the location of 
the zeros of A(^i,k). But this is immediately clear from (73). The integral has only 
poles at the zeros of A(^i , k). In case (I) there is one additional pole at Vi = v 2 +r| with 
residue (74a). The simple poles of G(vi ,v 2 ) at Xj +T|, where the Xj are the Bethe roots 
(see appendix A), are canceled by the simple poles of a(Vi,K) (see equation (A.6a)). 
In cases (II) and (IV) there is nothing to show. In case (III) we have one additional pole 
at Vi = v 2 — r| with residue (74b). The simple poles at Xj — r\ have vanishing residue 
due to (A.6a) and since 

res V| ^ G (v,,v 2 ) = -g^. (75, 



7 The exponential form - preliminary remarks 

The main result of [19] is the formula (1.12). It makes the calculation of arbitrary cor- 
relation functions possible, because the operators t* , b* , c* generate a basis of the space 
of quasi-local operators [5]. Although this formula proves the factorization of the cor- 
relation functions and allows, in principle, also for their direct numerical evaluation, 
it may be sometimes preferable to avoid the creation operators and to have an explicit 
formula for the correlation functions in the standard basis generated by the local op- 
erators . We believe that some form of the exponential formula discussed in the 
previous papers [2-4, 9] must be valid in case of temperature, disorder and magnetic 
fields as well. Unfortunately, the problem of constructing all operators that appear in 
this formula remains still open. We hope to come back to it in a future publication. 
Here we formulate the general properties we expect for these operators and show by 
examples how they should look like for short distances. 

From now on we shall use the notation and the terminology of the paper [9]. In 
particular we shall be dealing with the space W( a ) of quasi-local operators of the form 
q2as(0) q introduced there. First we define a density operator D* N : — > C which 
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generalizes that one defined by the formulae (33), (34) of [4], namely, for any quasi- 
local operator we define 

D* N (0) = (0} T .a, K (76) 

in such a way that 

D* N (e^...e m l?) =D N i\:±& u ...£ m \K,a) (77) 

where is the density matrix defined in (16). 
We expect that as before 

D* N (0) =tr a {exp(Q)(^ 2aS (°)0)} , (78) 

where tr a is the a-trace defined in [9] and where the operator £2 consists of two terms 
like that one constructed in [4] 

£2 = £2i+Q 2 - (79) 



In fact, the first term follows from [9, 19]. It must be of the form 



(80) 



^ 1= 7 2^C? / ^f|( ( °o(Ci/C2|a)-c (Ci,C2|K,cx))b(Ci)c(C2), 
where the function go was defined in [9], 

co (C|cx) = -(|^) 2 A^(C). (81) 



The second part in the right hand side of (79) should be of the form 

'2 



*2 = f^o mQm 



(82) 



where the operator t is yet to be determined. In some sense it must be the conjugate 
of the operator t*. The integration contour for both, Cli and Q.2, is taken around all 
simple poles £i , £2, £ = £j with j = 1, . . . ,m in anti-clockwise direction. The number 
m is the length of locality of the operator 0. 

Let us list some of the most important expected properties of the operator t. First, 
we expect that like t*(Q the operator t(Q is block diagonal, 

t(Q: Wa^Wo^, 

where, as was explained in [9], C is the space of quasi-local operators of 
spin s. We will deal below mostly with the sector s = 0. 

Then we expect t(Q to have simple poles at C, = Let us define 

t,- = res^.t(Q^- (83) 

while 

t*=t*(y. (84) 
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In contrast to (83) the operator t* is well defined only if it acts on the states X^q with 
/ < j. This will be always implied below. Let us denote /] (Q and respectively tj[k,i\ 
the operators defined on the interval [k, I] with k < j <l. 

We also expect that /^-matrix symmetry holds similar to the formula (2.16) of [9] 
fort*, 

s ; t M (Q =t M (Qs« for k<i<l. (85) 
Here as was defined in [9] 

s ; - = ^ + iI M+ i(^+i), (86a) 

%t+i&fe+i)( x ) = «M + ifeA i+ i)4y + i(^m)"' , (86b) 
R i j +1 (Q=P iJ+l R iJ+l (Q, (86c) 

where Kfj stands for the transposition of arguments and \j and P,j € End(V; <8> V/) 
is the transposition matrix. 
The further properties are: 

• Commutation relations 

[t,-,t t ]_ = [t J -,b(Ci)c(C2)]-=0. (87) 



• Projector property 

• Relations with t* 



t 2 j=tj. (88) 



tjt* k = t*t k for j^k, 

= t 5 t y = °- ( 89 ) 

• Reduction properties 

ti[i./](^%/])=^%/], (90a) 

t y[1)/] (^X [2) , ] ) = ^t 7 - [ 2 !/] (X [2!/] ) for 1< 7 </, (90b) 

t;[i,/](^[i,;-i]) = t y[ i,/_i](X [l!/ _ 1] ) for 1<J</, (90c) 

%,/](%/-!]) = 0. (90d) 

Let us comment on these relations. First of all the commutation relations (87) lead 
to the factorization of the exponential 

exp(£2) = exp(£2i +£^2) = exp(£2i)exp(£22) • (91) 

As we know (see [3, 4, 9] and earlier references therein) the operator £l\ becomes 
nilpotent when it acts on states of finite length. By way of contrast, the operator Q.2 is 
not nilpotent, but due to (87) and the projector property (88) one can conclude that 

exp(Q 2 )(<7 2aS( %, ;] ) = n(l -hhtl +P;ty[i ) /])(^[i,/])^ ( ° ) , (92) 
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where pj = p(^y). The reduction properties (90) look standard except for the first one. 
It is easy to see that we need all of them in order to have the reduction property (19b) 
of the density matrix, but still we do not have a good understanding of (90a). 

Let us show how the tj for s = explicitly look like in two particular cases, namely, 
for m = 1 and m = 2 where m = I — k + 1 and without loss of generality k = 1 . 
m = 1 : 

q° 



t 



■ a<fi (g> o\ 



(93) 



m 



2: 



tl[l,2] 



1 q a + q~ a 
4 q^-q-a 1 



1 



q-q 



-l 



q-q 



-1 



(a+G 2 



+ -q a M 



^a+1 -g-a-l 



•(g+g 2 -a^l) 



+ 



q a ~ l -q- a + l 



a 1 a+) 



+ 



12 



+ 



{q-q~ l ){q a + q- 



(? + ?- 1 )(9 a -?-°) 



2^ 



a+l 



-a-1 



-a+ 



T y(^ 1 /^-^i)-(aia 2 --ara 2 + ; 
j-^i A 2 + W^i) • Kc» 2 +ara+ 



(94) 



The operator 12(1,2] can be obtained from ti[i ,2] using the /^-matrix symmetry (85), 

-1 



t 



2[1,2] = Sltl[l,2]Si 



(95) 



One can check that all the above properties are fulfilled. 

It is interesting to understand how the limit a — > works, because, as one can see 
from (93)-(95), the operators ti and t2 are singular in this limit. More precisely, only 
the very first term in the expression (94) for tin^i is singular. In fact, they contribute 
into the density matrix only in such a combination that this singularity cancels. So, 
we actually need to calculate the residue with respect to a. In appendix B we will 
discuss this issue in more detail. We will also show how these operators are related to 
the 'fermionic' operators hj constructed in [2,4]. 



8 Conclusions 

The main result of this work is the description of the function go in terms of integrals 
involving the auxiliary function a and the function G, which are solutions of integral 
equations. As we have experienced in our previous work such type of description is 
efficient for performing the Trotter limit and for the actual numerical evaluation of 
correlation functions. We know from [19] that no other functions than go and p are 
required. Thus, together with [19], we have achieved a rather complete understanding 
of the mathematical structure of the static correlation functions of the XXZ chain. Our 
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results are equally valid in the finite temperature as in the finite length case, the only 
difference being a different driving term in the non-linear integral equation for the 
auxiliary functions o. 

We expect that our results open a way for further concrete studies of short range 
correlation functions at finite temperature as initiated in [2-4]. We hope that in the 
future it will also prove useful in studying field theoretical scaling limits as well as the 
large distance asymptotics of correlation functions in the XXZ chain. 

We have obtained our expression for go through a novel multiple integral represen- 
tation of the density matrix of the XXZ chain including a disorder parameter a. We 
think that this multiple integral representation is also interesting on its own right. 

We would like to point out that we obtained a remarkably beautiful and simple 
characterization of the function *P on the inhomogeneous finite lattice through equa- 
tion (67) and the residua (74). The function *P is important because it becomes the 
transcendental part of oo in the Trotter limit. 

We further performed a case study, looking for an operator t, adjoint to t* that 
allows us to write the density matrix in an exponential form even in the presence of 
a disorder parameter and a finite magnetic field. We obtained explicit expressions for 
t for m = 1,2, but so far could not find a general construction behind it. We hope to 
come back to this latter point in the future. 
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Appendix A: Derivation of the multiple integral representation 

A multiple integral representation for the density matrix in the non-twisted case a = 
was derived in [13]. Here we shall only indicate which modifications are necessary to 
include non-zero a and otherwise refer the reader to that work. 
First note that 

(K+«|r E E ;'(^,K)...r^fe B ,K)iK) _ (K|rj : fe B ,K)...^^,K)iK+«) 

(K+a|n? =1 f&,K)|K) (K|rT7 =1 ^-,K)|K+a) ' • 

because of the symmetry of the /^-matrix with respect to transposition. We may there- 
fore start our calculation with 

(K|r«'(^,K)...r p ";(U,K)iK + «) 

<K|n? =1 *&,K)|K+a) ' ( ' } 

which brings us closer to the notation of [13]. 

The left and right eigenvectors (k| and |k + oc) can be constructed by means of 
the algebraic Bethe ansatz. They are parameterized by two sets {k} = {kjYf^ and 
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{/u} = {/Lij}^ of Bethe roots, which are special solutions to the Bethe ansatz equations 



q- 2K d{l j ) ! lli sh(l j -X k + r\) 



(A.3) 

for 7=1,... ,N/2. By a(k) and d(k) we denoted here the vacuum expectation values 
of the diagonal elements of T(Q. For its K- twisted version we shall reserve the notation 



for the matrix elements. Then the eigenvectors (k| and |k + a) are 

(K\ = ({X}\ = (0\C(h)..-C(X N/2 ), 
\K + a) = \i M })=B(ft 1 )...B{ / t N/2 )\0), 



(A.4) 

(A.5a) 
(A.5b) 



where (0| and |0) are the left and right pseudo vacuum states. 

With the solutions {X} and {//} of the Bethe ansatz equations we associate the 
auxiliary functions 



a(A) = o(A,k) = -r— I I — r — r — , 

a[k) l = \sh(X-X k -r[) 

a a (A) = a(A,K + a) = — I I — 



(A.6a) 
(A.6b) 



and the ratio of ^-functions 



|=i sh(X-^-) 

Then ({k}\ is the 'dominant' left eigenvector of t(C,,K) with eigenvalue 



(A.7) 



(l + o(A,)), 



(A.8) 



and similarly |{ju}) is the dominant right eigenvector of the oc-twisted transfer matrix 
t(C,,K + a) with eigenvalue 



A(^,K + oc) =q K+a a(X) 



|/ 2 sh(X-^-T!) 



.|=1 sh(X-/ly) 

Dividing (A.9) by (A.8) we obtain the identities 



(l + a„(X)) (A.9) 



P(C) 



l + o(X) 



. q a W-i\)r 1 w = 



l + a(X) 



^XX + ri)^- 1 ^) (A. 10) 
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which will be needed below. Here we have introduced the notation o = 1 /a and a a = 
l/a a . In the derivation of the multiple integral formula below we shall use that p is 
analytic and non-zero inside the canonical contour C which follows from the Bethe 
equations and from the explicit form of the vacuum expectation values a and d. 

The derivation of the density matrix for a = in [13] is divided into two steps. 
Step 1 is the derivation of the 'general left action' ({X}\T^ [ (^,k) . . . r p a "(^ m ,K) of 
a string of monodromy matrix elements on the left dominant eigenvector. This step 
remains the same as before. The general left action is given by Lemma 1 of [13]. In a 
second step one must calculate ratios of scalar products of the form 



({v+}u{X-}|M) 

(WlM)n7=iA(^,K) 



(A. 11) 



which are generated in step 1. A useful expression for these ratios in the untwisted 
case a = is provided by lemma 2 of [13]. This needs to be modified here. 

The notation in (A. 11) is meant as follows. We divide the sets {X} and {v} = 
{ y j}"j=\ = {^j}J=i mto disjoint subsets {X + }, {X~} and {v + }, {v~}, such that their 
unions are {X + } U {X~} = {X} and {v + } U {v~} = {v}. The number of elements in a 
set {x} will be denoted \x\. We shall assume that \X + \ = [v + 1 = n. For a given partition 
of {X} into {X + }, {X~} we order the Xs such that 



(Xi,...,X N / 2 ) - {X 1 ,...,X n ,...,X N , 2 _ n 



(A. 12) 



and we define 



(Xi,...,X N / 2 ) - (Vj ,...,V„ ^ ,...,X N/2 _ n ). 



(A.13) 



Then using lemma 2 of [14] (with the roles of X and (i interchanged) we arrive after 
some trivial cancellations at 



y=Wv y )(l + a(v ; .)) 



]X+| 

n 

j= i 



1,^X1 



b(x k -xjy 



a(Xj)(l + a(vj)) 



det 



X+\ 



c(Xj 



iHXj 

det N/2 N(fij,X k ) lj=i 



(A. 14) 



where 



N(nj,X k ) = e(jij - X k ) - e(X k - /jj)a a (X k ) . 



(A. 15) 



In the following we concentrate on the term E in the second line of (A. 14). We 
want to transform it into a form that allows us to perform the Trotter limit jV — > °°. We 
define column vectors k = 1, . . . , |v + | = n, and k = 1, . . . ,N/2, by 

(u k ) j = 4>(v+ -Ti)% y ,v+) = *(v+ -Ti)(e(/i 7 - -v+) - e(v+ -Hj)Oa(V+)) , 



(y k ) J =q>(X k -r\)N(nj,X k ) =e(nj-X k )<b(X k -r\) + q 2a e(X k -Hj)<b(X k + T]) . 



(A. 16) 
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In (A. 16) we used the Bethe ansatz equations (A. 3) to eliminate q 2K d(X k )/a(X k ). We 
have 

„ _ det(m, . . . ,u„, v„ + i , . . . , 



det(\ u ...,\ N/2 ) 



(A.17) 



Next we use a trick we learned from N. Kitanine [20] and which in a similar form 
originally appeared in [16]. Define a matrix X with matrix elements 

Xj = q a cth^ l -X j )re^ i r\X) (A.18) 

and column vectors U k = U(v^) = Xu k , and = X\ k . Then 

„ _ det(Ui,...,U n ,V„ + i,...,V^/ 2 ) 



detrVi,...^) 



(A. 19) 



The vectors U k and \ k are easily calculated by means of the residue theorem. Choose 
two simple closed contours and 6^ such that all Bethe roots Xj are inside Sx but 
outside and all Bethe roots \ij are inside but outside S^. Then 

(v*y = 

N/2 

£ 4 a cth( w - Xj) res x=w <\>~ l (X) [e( W - X k )<\>(X k + q~ 2a e(X k - w)<\>(h + 
i=i 

= [ ^cthCu-X y )r 1 ()")^ <x eCu-^(^-Ti) + ? - a e(^- / i)^ + Ti)] . 

J Gfj £iw s ^, ' 

= M 

(A.20) 

The function / is periodic with period i7t, and limR e// ->±<»/(/i) = 0. Moreover, (p) 
is analytic inside Qx- Hence, for j / k, 



(VicY = ~J e Iri-fM = ~( re& v=h + Tes n=h + res »=h+r\ +res /J=Xjt _ T1 )/(//) 



: q- a cth(Xj -X k -r\)- q a cih{Xj -X k + r\) = K a (Xj - X k ) . (A.21) 



For j = & we have a non-trivial residue from the second order pole at Xj, and 

= (q-^(Xj + T]) - q^Xj - Tl))3^- 1 (m) \ p=Xj + K a (0) . (A.22) 
Similarly, we obtain 

(U k y = q<*cth(Xj - v+)r 1 «;H>(v+ -n) [1 + a„(v+)] - q"cth(Xj - v+ +T1) 

-</*cth(^ - v+ -ri)r !(v+ +Tl)*(v+ -Ti)a a (v+) . (A.23) 

Following once more [20] we eliminate the function (]) from (A.22) and (A.23). 
This can be done by means of the identities (A. 10). For the first term on the right hand 
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side of (A.22) we obtain 



lim fl '(X)p(C) 

A — 



a(X) 



1 



l + a a (X) l + Oa(X) 



-o'(^-)p(Cy), (A-24) 



where we used (A. 10) in the second equation and the Bethe ansatz equation a(Xj) = — 1 
in the third equation. It follows that 



Eliminating (j) from (A.23) by means of (A. 10) we end up with 



(A.25) 



(U,y=cth(^-v+)(l + a(v+))p(^) 

- <7 a cth(^ - v+ + r\) - ^ a cth(X y - v+ - ri)a(v+) . (A.26) 

From here on we can proceed as in [14]. Define a matrix V = (\\,...,\ N / 2 ) an d 
a vector W(v) = V _1 U(v). Then 

S = det N/2 (V _1 Ui , . . . , V~ l \] n ,e n+ i ,e N / 2 ) = 

det„((e 7 -,y- 1 U,}) = det„(W(v+y) . (A.27) 

W(v^) is the solution of the linear equation VW(Vj) = U(v^), or, explicitly, 

cth(^-v+)(l + a(v+))p(^)-^cth(^-v+ + Ti)-, ? -«cth(^-v+-Ti)a(v+) 

N/2 

= -W(v+)V(^-)p(g) + £ K a (Xj - X/)W(v+)< . (A.28) 
This can be transformed into a linear integral equation. For this purpose define 

q a C\h{X-V + l\) q~ a c\h{X - V - X\) 
Cr(A,,VJ = — — \ 



1 + o(v) 



This function is defined such that 



l + c(v) 

c^-v) P@+ "£ ^,: x l w(v> ' . (a.29, 

i=i 



1 + o(v) 



G(Xj,v+) 



pfo)a'(A,j)W(v+y 
l + a(v+) 



(A.30) 



We shall assume that v is located inside the canonical contour C shown in figure 1. By 
construction G(X,v) is then meromorphic inside C and has a single simple pole with 
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residue — p(^) at X = V. Using (A.30) it follows that 
G(X, V+) = q- a Cth(l - V+ - T!) - Cth(X - V+)p(£) 



+ 



e 2ni p(e")(l + cOi)) 



(A.31) 



which is a linear integral equation for G. 

Combining (A.27) and (A.30) we infer that 



ffl i + q(y 7 + ) 

LH tf(*?)p(ff) 
Inserting this into (A. 14) we arrive at 



det| X+ |(G(A,+ v+)). 



(A.32) 



y«(v7)(i+o(v7)) 



n 



det| X+ |(G(A,+ v+)) 



det 



|X+| 



*(V- V *) 



(A.33) 

This replaces the expression in lemma 2 of [13]. Comparing the two expressions we 
see only one explicit difference which is the appearance of the additional factor p(£j~) 
in the denominator. It always comes with a factor a'(Xj~). Therefore it is easy to trace 
the modification required in the derivation of the multiple integral formula in [13]. 
There is, however, also an implicit change in the formula. The residue of G(k,v) at 
X = V is — p(^) instead of —1. This causes another tiny modification in the derivation 
of the multiple integral formula. In the first equation (65) of [13] additional m — 
factors of p _1 appear which together with the explicit factors gives exactly one 
factor per integral. The final result for the density matrix matrix is 



(K\T^g u K)...T^(t >m ,K)\K + a) 
<K|IT? =1 ^;,K)|K + a) 

(| ."/ F M) 



fr f FM fr f cfrj PjM 

}J Je 2ni p(e^)(l + a(^-)) J L . =| y +1 Je 2ni p(e^)(l + afrj 

dfSt[-G{fij,\ k )] 



)) 



sh( A / i - w -ri)sh(v (t -Vy) 

!<}<k<m 



(A.34) 



where we have used the notation of [13], i.e. 



Fj(X)= nsh(X-v*-Ti) J] sh(A,-v*), j=l,...,\a+\, (A.35a) 

k=\ k=xj+l 

Fj(X)= nsh(X-v, + ri) J] sh(X-Vfc), j = \a + \ + 1,. . . ,m, (A.35b) 

k=l k=xj+l 



28 



and for j = 1, .. . , |oc + | we define xj to be the position of the (|oc + | — j+ l)th '+' in the 
sequence of upper indices ((X/)>=i while for j = \a + \ + 1, . . . ,m it means the position 
of the (j- |oc + |)th '-' in the sequence of lower indices (Py)^^ 

Finally we replace by ^ m - k +i and define = a m - k+ i, e[ = $ m -k+i for = 
1, . . . ,m. We denote the position of the jth '+' in (£/)yLi by £+ , the position of the 
jth '-' in (e^)7=i by ej. Then 

, y =( m - e ' + + 1 , i=1 '-;" , (a.36) 

L m_e m-j+i + 1 J=P + h---,m 

where p = |oc + | 

i = 1 ''"'/ , (A-37) 
j =p+l,...,m 

we obtain Xj = m — £j + 1 , j = 1 , . . . , m. Finally, setting (p) = Fj (p) and Ff (p) = 
F and using (A.l), we arrive at (29). 




Appendix B: Relation with previous results 



In our previous work [2,4], before we knew the multiple integral representation for 
finite a, we conjectured formulae which we claimed to hold in the limit a — > 0, relevant 
for the physical correlation functions of the XXZ chain. We introduced a function, 
say, C0 id(jui,/i2;cc) defined by an integral formula involving, among other functions, 
a function G id which is different from G defined in (31). Here we explain why our 
previous results remain unaltered in the limit a — > if we replace (floidG"ijj"2;oc) with 
— 0)(^i,^2|K,a) + COo(^|a) (the minus sign is due to a change of conventions in which 
we followed [9, 19]). 

Our previous ad hoc definitions were 

C0oid(vi,v 2 ;oc)-coo(^a) = -^ ax F i d (v 2 ,Vi;-a)+A\|/(^), (B.l) 

^oid (v 2 , v i ; -a) = 2 J dm Q {p) G oXA (p, v 2 ; - a) (q a cth {p - v i - x\ ) - cthfoi - v i ) ) , 

where dmoip) = dm(p) | a= o and where G id was the solution of the linear integral equa- 
tion 



G id(X,v;-a) = 

^ a cth(X-v-ri)-cth(X-v)+ f dm (n)K a (X-v)G{n,v;-a). (B.2) 

Je 

Comparing with (52), (50) and (31) we see that apart from some conventional 
sign changes the only difference is that the function p is replaced by unity in the old 
definitions. Using that (0o(^|0) = we conclude in particular that 



a>(£i,$2|K,0) = -co id(vi,v 2 ;0) 



(B.3) 
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and that G(k,v)\ a =o = G o id(^,v;0) = Go(k,v) which satisfies the integral equation 

G (X,v) =e(v-A,)+ f dm (M)K(X-v)G (M,v). (B.4) 
This integral equation implies the symmetry 

©($l,$2|K,0)=<0($2,$i|K,0). (B.5) 

Let us define 

co'(^^ 2 ) =3 a (r a co^i^2|K,a))| a=0 , (B.6a) 

»dd( v i» v 2) =9a(^ H cOoid(vi,V2;a))| a=() . (B.6b) 

We shall show below that 

<u(Vi,v 2 ) = -co'(^,^)) . (B.7) 

Accepting this for a moment let us insert -co + COo instead of oo id into our previous 
formula for the exponential form, e.g. into equation (35) of [2], which actually means 
to use (80) for the t independent part of the exponential form. Then using (B.3), (B.5) 
and (B.7) and the fact that ^ _a COo(^|a) = 0(oc 2 ) we recover our previous result, namely 
equation (37) of [2], in the limit a — > 0. 

It remains to prove (B.7). For this purpose consider 

I((Oold(Vl,V2)+(0 / (^1^2)) = (P'GO-P'fe))^) 

-p'(£i) / dmo(/i)Go(/i,V2)cth(/i-Vi)- / dm Gu)pV)GoGu,v 2 )e(vi - M ) 
Je Je 

+ f dm (A()G / (A ( ,V2)e(v 1 -^), (B.8) 

where p'(£) = 3 a p(^)|o=o, and G'(/i,v 2 ) = a a (G( i u,v 2 ) + G old (X,v;a))| a=0 . Taking 
the a-derivative of (31) and (B.2) we find 

G / (X,v 2 ) = -p / fe)cth(X-v 2 )- /" dm ( Af )p / (e")Go(A/,V2)^- A( ) 

Je 

- /" dmoCu^^-^G'^Vz). (B.9) 

J G 

Using (B.4) and (B.9) we can eliminate G' from (B.8) by means of the 'dressed func- 
tion trick' . We arrive at 

K«old(Vl,V 2 )+(0 / (^^2)) = (p'^O-p'fe))^) 

-p'(Si) / dmoCu)Go(/i,v 2 )cth0u-Vi)-p / (^2) / dm Cu)G (/i,Vi)cth(/i- v 2 ) 
Je Je 

- f dmoCu)p'(e^)Go(/i,v 1 )GoCu,V2) (B.10) 
Je 
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which is obviously symmetric. Then (B.7) follows if one takes into account that 
ro oid(^i'^ 2 ) * s antisymmetric, which was shown in [4]. 

Now let us discuss how the operators tj described in the section 7 are related to 
the operators hj (see formulae (40)-(42) of [2]) in the limit a — > 0. As we mentioned 
above the operators tj have a pole of first order when a — > 0. On the other hand, it 
follows from the formula (92) that the density matrix depends only on the combination 
(Pj — l)tj. Since lim a ^oP; = 1, 1 — p 7 = 0(a), and one obtains for (pj — l)ty a finite 
result in the limit a — > 0. So we actually need only the residues of tj. Let us define 



t 



(0) 



lim(l 

cc^O 



(B.ll) 



Since for the moment we do not have an explicit formula for tj in the general case, let 
us describe the relation with hj again only for the cases m = 1 , 2 where we know the 
explicit result (93) and (94), (95). It is enough to consider j = 1. 
For m = 1 



t 



(o) 

1[U] 



-h 1[M] = -^/c 



For m = 2 we simply get from (94) 



(B.12) 



*(0) 

l l[l,2] 



1 

-4" 



q-q 



(<J+G 2 -^0+) 



(B.13) 



and 



.(0) . , 

l l[l,2] + h l[1.2] 



o\o z 2 



'+( 



^2 + ^1 



(a+a 2 



(B.14) 



It is interesting to note the following. In spite of the fact that the operators hj 
are fermionic and the tj are bosonic, the limit of the oc-trace of the corresponding 
exponentials should coincide 



lim tr a {exp(£logp;t;) 



exp(-£(p 7 h,)}(^ s (°)0)=0 (B.15) 



where cpy = (p(^-|K,0). We do not have a complete understanding of this relation for 
the moment, but we can follow how it works for m = 1,2. The case m = 1 is trivial 
while the case m = 2 is more interesting. For instance, one could expect that terms 
containing cpiq>2 appear, but they do not. For the fermionic formula this is a simple 
consequence of the anti-commutation relations of the operators hj. The reason why 
do they not appear for the bosonic formula is different, namely, because of the relation 



t (0) t (0) 
future. 



0, which can be checked. We plan to study this point more carefully in the 
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